#ifndef PBVP_VOLUME3D_CRITICAL_H
#define PBVP_VOLUME3D_CRITICAL_H

#include "../data/reg3dvolumedata.h"
#include "tetra.h"
#include "contreeconstructor3d.h"
#include "../mact/moment.h"

namespace PBVP {
	/// A volumetric data class for constructing CT and DCT
	/**
	 * This class divide each regular cell of 3D volumetric data
	 * into six tetrahedra
	 */
	class Volume3DCritical {
	public:
		/// Constructor
		Volume3DCritical(Reg3DVolumeData* pData);

		///
		virtual ~Volume3DCritical() {}

		/// Get the values and coordinates of the four vertices of the tetrahedron of tid
		/**
		* 
		* @param tid: id of the tetrahedron
		* @param vals: return values
		* @param x0: return coordinates of vertex 1
		* @param x1: return coordinates of vertex 2
		* @param x2: return coordinates of vertex 3
		* @param x3: return coordinates of vertex 4
		*/
		inline void getTetraData(int tid, float vals[4], float x0[3], 
								 float x1[3], float x2[3], float x3[3]);


		/// Get the ids of the 4 vertices of a tetrahedron
		/**
		 * @param tid Id of the tetrahedron
		 * @param vids Returned vetex Ids
		 */
		inline void getTetraVertIds(int tid, int vids[4]);

		/**
		*	Return a tetrahedron that contains the given edge
		*	@parameter 
		*/
		int getTetraIndex(int idx1[3], int idx2[3]);

		/// return the np-th property function value at index (i, j, k)
		/**
		* 
		* @param i
		* @param j
		* @param k
		* @param np The np-th property. The property index starts at 1
		* 
		* @return float
		*/
		float getPropertyValue(unsigned int i, unsigned int j, unsigned int k, int np = 1) const {
			assert(np <= getNumofProperties());
			return m_pData->getValue(i, j, k, np); 
		}

		/// return the np-th property function value at coord[3]
		/**
		* 
		* @param coord
		* @param np The np-th property. The property index starts at 1
		* 
		* @return float
		*/
		float getPropertyValue(float coord[3], int np = 1) const {
			assert(np <= getNumofProperties());
			return m_pData->getValue(coord, np);
		}

		/// get the number of properties functions
		/**
		* The first variable in the volumetric data is the shape
		* function. The rest are property functions. 
		*/
		int getNumofProperties() const {
			return m_pData->numOfVariables() - 1;
		}

		/**
		* Check if two vertices are adjacent to each other
		* @param v1
		* @param v2
		* 
		* @return bool
		*/
		inline bool areNeighbors(int v1, int v2);

		/**
		* Check if two vertices are adjacent to each other
		*
		* @param idx1
		* @param idx2
		* 
		* @return bool
		*/
		inline bool areNeighbors(int idx1[3], int idx2[3]);

		/**
		* get the neighbors of a tetrahedron
		* @param side: 
		* @return the neighboring tetrahedron id 
		*/
		int tetNeighbor(int tid, int side)
		{
			int cid = tid / 6;
			int tet = tid % 6;
			int idx[3];

			m_pData->cell2index(cid, idx);

			idx[0] += Tetra::six_tetra_neighbors[tet][side][0];
			idx[1] += Tetra::six_tetra_neighbors[tet][side][1];
			idx[2] += Tetra::six_tetra_neighbors[tet][side][2];

			if (!isValidCell(idx)) return -1;
			return(6*m_pData->index2cell(idx[0], idx[1], idx[2]) + Tetra::six_tetra_neighbors[tet][side][3]);
		}

		/**
		* Find a tetrahedron edge that intersects a cut value 
		*  
		* @param cut_val
		* @param range
		* @param v1 Start vertex to search from
		* 
		* @return int
		*/
		int getCutEdge(float cut_val, float range[2], int &v1);

		int getCutEdge(float cut_val, int &v1);

		/// Get the n-th property values on four vertices of a tetrahedron
		void getTetraProperty(int tid, float pot[4], int np = 1);

		/** 
		*  Check if two ranges intersect each other.	
		*/
		bool areDisjointRanges(float rang1[2], float rang2[2]) 
		{
			if (rang2[0] > rang1[1]) return true;
			if (rang2[1] < rang1[0]) return true;
			if (rang1[0] > rang2[1]) return true;
			if (rang1[1] < rang2[0]) return true;
			return false;
		}
	private:
		Reg3DVolumeData *m_pData;
		///
		unsigned int m_Dim[3];
		///
		float m_Span[3], m_Orig[3];

		/// Check if the index is valid cell index
		bool isValidCell(int idx[3])
		{
			return(idx[0] >= 0 && idx[1] >= 0 && idx[2] >= 0 &&
				idx[0] < (m_Dim[0]-1) && idx[1] < (m_Dim[1]-1) && idx[2] < (m_Dim[2]-1));
		}

		/// Check if the index is valid vertex index
		bool isValidVertex(int idx[3])
		{
			return(idx[0] >= 0 && idx[1] >= 0 && idx[2] >= 0 &&
				idx[0] < m_Dim[0] && idx[1] < m_Dim[1] && idx[2] < m_Dim[2]);
		}

		/**
		* Sort tetra vertices by their values and remove any degeneracy 
		*/
		void sortVerts(float* &x0, float* &x1, float* &x2, float* &x3, float v[4], float f[4]);

		/**
		* Find the neighboring vertex of the maximum value
		* @param vid
		* 
		* @return int
		*/
		int getMaxValNeighbor(int vid);

		/*
		*	Find the neighboring vertex of the minimum value
		*/
		int getMinValNeighbor(int vid);


		float getCPValue(int id) const {
			return m_pData->getValue(id, 0);
		}

		/**
		*  Calculate the volume of a connected component with given range
		*  @param fint Integral of a second function on the volume
		*/
		float calcVolume(float range[2], int start_tid, float &fint);

		/**
		* Calculate volumetric moments for a shell of given range
		*/
		VolMoments calcMoments(float range[2], int start_id);

		/*
		*	The volume of the reg3D data
		*/
		float getVolume() const {
			// assumed cube cells
			return(m_Dim[0]-1)*(m_Dim[1]-1)*(m_Dim[2]-1)*m_Span[0]*m_Span[1]*m_Span[2];
		}

		/**
		* Calculate Carbo index for a given orientation
		*/
		//friend float CarboIndex(const ConTreeConstructorReg3D& reg1, float **R1, float *C1,
		//	const ConTreeConstructorReg3D& reg2, float **R2, float *C2);

		///
		float subFuncInteg(float x1[3], float x2[3], float x3[3], float x4[3], float v[4], float f[4], float range[2])
		{
			float *p1 = x1, *p2 = x2, *p3 = x3, *p4 = x4;
			float vals[4], func[4];
			for (int i = 0; i < 4; i++)
			{ // make the copies so we don't mess with the original 
				vals[i] = v[i];
				func[i] = f[i];
			}
			sortVerts(p1, p2, p3, p4, vals, func);
			Tetra tetra(p1, p2, p3, p4, vals[0], vals[1], vals[2], vals[3]);
			return(float)(tetra.funcIntegral(func[0], func[1], func[2], func[3], range[1]) -
				tetra.funcIntegral(func[0], func[1], func[2], func[3], range[0]));
			//return tetraFuncIntegral(p1, p2, p3, p4, vals[0], vals[1], vals[2], vals[3], func[0], func[1], func[2], func[3], range[1]) - 
			//	   tetraFuncIntegral(p1, p2, p3, p4, vals[0], vals[1], vals[2], vals[3], func[0], func[1], func[2], func[3], range[0]);
		}

		void subFuncIntegMul(float x1[3], float x2[3], float x3[3], float x4[3], float v[4],
			float *f1, float *f2, float *f3, float *f4, 
			float *out, int n, float range[2])
		{
			float *ou = new float[n];
			float *ol = new float[n];
			/*float *p1 = x1, *p2 = x2, *p3 = x3, *p4 = x4;
			float vals[4], func[4];
			for(int i = 0; i < 4; i++) { // make the copies so we don't mess with the original 
			vals[i] = v[i];
			}
			sortVerts(p1, p2, p3, p4, vals, func);
			float *g1 = f1, *g2 = f2, *g3 = f3, *g4 = f4;
			float vals2[4], func2[4];
			for(int i = 0; i < 4; i++) { // make the copies so we don't mess with the original 
			vals2[i] = v[i];
			}
			sortVerts(g1, g2, g3, g4, vals2, func2);

			tetIntegMultiFuncs(p1, p2, p3, p4, vals[0], vals[1], vals[2], vals[3], g1, g2, g3, g4, ou, n, range[1]);
			tetIntegMultiFuncs(p1, p2, p3, p4, vals[0], vals[1], vals[2], vals[3], g1, g2, g3, g4, ol, n, range[0]);
			if(!(out[0] <= 1)) {
			printf("---vals = [%f %f %f %f], range = [%f %f], vol = %f, ou = %f, ol = %f\n", 
			vals[0], vals[1],
			vals[2], vals[3], range[0], range[1], out[0], ou[0], ol[0]);
			printf("p1 = (%f %f %f)\n", p1[0], p1[1], p1[2]);
			printf("p2 = (%f %f %f)\n", p2[0], p2[1], p2[2]);
			printf("p3 = (%f %f %f)\n", p3[0], p3[1], p3[2]);
			printf("p4 = (%f %f %f)\n", p4[0], p4[1], p4[2]);
			}
			*/
			Tetra tetra(x1, x2, x3, x4, v[0], v[1], v[2], v[3]);
			tetra.multiFuncsIntegral(f1, f2, f3, f4, ou, n, range[1]);
			tetra.multiFuncsIntegral(f1, f2, f3, f4, ol, n, range[0]);
			for (int i = 0; i < n; i++)
			{
				out[i] = ou[i] -ol[i];
				//printf("%d %f\n", i, out[i]);
			}
			//assert(out[0] >= 0);
			delete[] ou;
			delete[] ol;
			return;
		}   

		float subVolume(float x1[3], float x2[3], float x3[3], float x4[3], float v[4], float range[2])
		{
			float *p1 = x1, *p2 = x2, *p3 = x3, *p4 = x4;
			float vals[4], func[4];
			for (int i = 0; i < 4; i++)
			{ // make the copies so we don't mess with the original 
				vals[i] = v[i];
			}
			sortVerts(p1, p2, p3, p4, vals, func);
			Tetra tetra(p1, p2, p3, p4, vals[0], vals[1], vals[2], vals[3]);
			return(float)(tetra.innerVolume(range[1]) - tetra.innerVolume(range[0]));
			//return tetVolume(p1, p2, p3, p4, vals[0], vals[1], vals[2], vals[3], range[1]) - 
			//	   tetVolume(p1, p2, p3, p4, vals[0], vals[1], vals[2], vals[3], range[0]);
		}
	};

	void Volume3DCritical::getTetraData(int tid, float vals[4], float x0[3], 
										float x1[3], float x2[3], float x3[3])
	{
		// cell id
		int cid = tid / 6;
		int tet = tid %6;
		int idx[3];
		m_pData->cell2index(cid, idx);

		vals[0] = m_pData->getValue(idx[0]+Tetra::six_tetra_index[tet][0][0], 
									idx[1]+Tetra::six_tetra_index[tet][0][1], idx[2]+Tetra::six_tetra_index[tet][0][2]);
		vals[1] = m_pData->getValue(idx[0]+Tetra::six_tetra_index[tet][1][0], 
									idx[1]+Tetra::six_tetra_index[tet][1][1], idx[2]+Tetra::six_tetra_index[tet][1][2]);
		vals[2] = m_pData->getValue(idx[0]+Tetra::six_tetra_index[tet][2][0], 
									idx[1]+Tetra::six_tetra_index[tet][2][1], idx[2]+Tetra::six_tetra_index[tet][2][2]);
		vals[3] = m_pData->getValue(idx[0]+Tetra::six_tetra_index[tet][3][0], 
									idx[1]+Tetra::six_tetra_index[tet][3][1], idx[2]+Tetra::six_tetra_index[tet][3][2]);

		for (int i = 0; i < 3; i++) {
			x0[i] = m_Orig[i]+(idx[i] + Tetra::six_tetra_index[tet][0][i])*m_Span[i];
			x1[i] = m_Orig[i]+(idx[i] + Tetra::six_tetra_index[tet][1][i])*m_Span[i];
			x2[i] = m_Orig[i]+(idx[i] + Tetra::six_tetra_index[tet][2][i])*m_Span[i];
			x3[i] = m_Orig[i]+(idx[i] + Tetra::six_tetra_index[tet][3][i])*m_Span[i];
		}
	}

	void Volume3DCritical::getTetraVertIds(int tid, int vids[4])
	{
		// cell id
		int cid = tid / 6;
		int tet = tid %6;
		int idx[3];
		m_pData->cell2index(cid, idx);

		vids[0] = m_pData->index2vert(idx[0]+Tetra::six_tetra_index[tet][0][0], 
									  idx[1]+Tetra::six_tetra_index[tet][0][1], 
									  idx[2]+Tetra::six_tetra_index[tet][0][2]);
		vids[1] = m_pData->index2vert(idx[0]+Tetra::six_tetra_index[tet][1][0], 
									  idx[1]+Tetra::six_tetra_index[tet][1][1], 
									  idx[2]+Tetra::six_tetra_index[tet][1][2]);
		vids[2] = m_pData->index2vert(idx[0]+Tetra::six_tetra_index[tet][2][0], 
									  idx[1]+Tetra::six_tetra_index[tet][2][1], 
									  idx[2]+Tetra::six_tetra_index[tet][2][2]);
		vids[3] = m_pData->index2vert(idx[0]+Tetra::six_tetra_index[tet][3][0], 
									  idx[1]+Tetra::six_tetra_index[tet][3][1], 
									  idx[2]+Tetra::six_tetra_index[tet][3][2]);
	}

	bool Volume3DCritical::areNeighbors(int v1, int v2)
	{
		int idx1[3], idx2[3];
		m_pData->vert2index(v1, idx1);
		m_pData->vert2index(v2, idx2);

		for (int i = 0; i < 14; i++)
		{
			if ((idx1[0]+Tetra::six_simplex_neighbors[i][0] == idx2[0]) &&
				(idx1[1]+Tetra::six_simplex_neighbors[i][1] == idx2[1]) &&
				(idx1[2]+Tetra::six_simplex_neighbors[i][2] == idx2[2]))
				return true;
		}
		return false;
	}

	bool Volume3DCritical::areNeighbors(int idx1[3], int idx2[3])
	{
		for (int i = 0; i < 14; i++)
		{
			if ((idx1[0]+Tetra::six_simplex_neighbors[i][0] == idx2[0]) &&
				(idx1[1]+Tetra::six_simplex_neighbors[i][1] == idx2[1]) &&
				(idx1[2]+Tetra::six_simplex_neighbors[i][2] == idx2[2]))
				return true;
		}
		return false;
	}
}
#endif 